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ABSTRACT 

Scale-free disks have no preferred length or time scale. The question has been raised 
whether such disks have a continuum of unstable linear modes or perhaps no unstable 
modes at all. We resolve this paradox by analysing the particular case of a gaseous, 
isentropic disk with a completely flat rotation curve (the Mestel disk) exactly. The 
heart of the matter is this: what are the correct boundary conditions to impose at the 
origin or central cusp? We argue that the linear stability problem is ill-posed. From any 
finite radius, waves reach the origin after finite time but with logarithmically divergent 
phase. Instabilities exist, but their pattern speeds depend upon an undetermined phase 
with which waves arc reflected from the origin. For any definite choice of this phase, 
there is an infinite but discrete set of growing modes. Similar ambiguities may afflict 
general disk models with power-law central cusps. The ratio of growth rate to pattern 
speed, however, is independent of the central phase. This ratio is derived in closed form 
for non self-gravitating normal modes. The ratio for self-gravitating normal modes is 
found numerically by solving recurrence relations in Mcllin-transform space. 

Key words: instabilities - hydrodynamics - galaxies: kinematics and dynamics - 
galaxies: spiral 



1 INTRODUCTION 



A scale- free disk, whether gaseous or stellar, has no charac- 
teristic frequency. It is therefore puzzling how such a disk 
can distinguish a characteristic frequency for its marginally 
stable modes. This paradox was first pointed out by Zang 
(1976) in a famous Massachusetts Institute of Technology 
Ph. D. thesis (supervised by Alar Toomre), which explored 
the stability of the scale-free stellar dynamical disks with 
completely flat rotation curves. Zang showed that if a scale- 
free disk admits one unstable mode, then a two-dimensional 
continuum can be constructed by self-similar scaling. Think- 
ing this improbable, Zang argued that it is more likely that 
scale-free disks admit no normal modes whatsoever. If true, 
this is astonishing, as it must hold good even for the com- 
pletely cold disk. Subsequently, the problem was re-visited 
by Lynden-Bell & Lemos (1993), who took the opposite view 
and argued that all the modes become unstable together. 
The possibility that an entire continuum of modes, each with 
a different frequency, might all become simultaneously un- 
stable had earlier been envisaged by Birkhoff in his classic 
1960 textbook on Hydrodynamics. If true, this also is aston- 
ishing, as the disks would become unstable at all scales to 
normal modes of all frequencies. 

This paper resolves the paradox by showing how to 
construct an infinite but discrete set of normal modes. Sec- 
tion 2 introduces the disks under scrutiny and sets up the 



mode equations for the perturbed quantities. Section 3 ex- 
amines the easier but still ambiguous problem of modes 
without self-gravity, whereas Section 4 deals with the fully 
self-gravitating case. 



2 EQUILIBRIA AND PERTURBATIONS 
2.1 Equilibria 

Our disks are scale-free in every respect. For definiteness, we 
suppose that they are razor-thin and fully self-gravitating in 
equilibrium. In fact, these properties are inessential as long 
as there is no preferred scale. Wherever possible, we use 
the same notation as in the recent stability analysis of the 
stellar dynamical power-law disks carried out by Evans & 
Read (1998a,b). Some basic formulae are repeated here for 
convenience and completeness. 

In equilibrium, the surface density and self-consistent 
potential are 



Scq(-R) — £ 



y>c q (i?) 



lid 

v l (1LY S 

8 KRoJ ' 



(1) 



in which vg is a constant reference velocity: 
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' r[i(i + /3)]r [|(2-/3)] 

Here, t/j is a positive quantity; the gravitational acceleration 
is V?/). As (3 — *■ 0, w| — > 27rG'So_ Ro = tip, and the potential 
becomes that of the Mestel disk ( Mestel 1963j ) , 

(*) ■ 



The reference radius i?o and surface density Eo appear for 
the sake of dimensional consistency but do not introduce 
preferred scales, since the replacements 



Ro -> sR , E -> s E , u fl -> s ^Sg, (4) 

preserve the form of equations [except in the case = 

0, where the potential (^) gains a physically unimportant 
additive constant]. 

The disks studied here are gaseous and isentropic, un- 
like their stellar dynamical counterparts explored by Evans 
& Read ( 1998a, b). The axisymmetric and the neutral modes 
of the gaseous power-law disks have already been the sub- 
ject of recent attention (e.g., Schmitz & Ebert 1987; Lemos, 
Kalnajs & Lynden-Bell 1991; Syer & Tremaine 1996). The 
square of the sound speed varies with radius as 



1 (no) 



(5) 



For convenience, let us introduce a dimensionless "temper- 
ature" 



e : 



\d%l) e Jd\nR\ 



(6) 



Radial hydrostatic equilibrium requires the square of the 
circular velocity to be 



*4c = (i-e)^(|-) . 

The Mach number is defined as: 
M 



C 



i - e 
e ■ 



The angular velocity and epicyclic frequency, 

1/2 



Q : 



R 



1 d tu ^2 
RSdR {Rv - ] 



(T) 



(8) 



(9) 



(10) 



both scale as R 1 , and their ratio is 

We assume — 1 < (3 < 1 so that k is always real. Toomre's 
con ditio n Q > 1 for local axisymmetric stability (Toomrc 
19640 related to 9 by 



Q 



4 



^(2-/3)9(1-6). 



ttGE V27rGE -R 0/ 
In terms of the Mach number, this becomes: 
( vl \ M 



= 2^2~Ji 



\2nGT, R J 1 +M 



(11) 



(12) 



The factor in parentheses, which can be read off from equa- 
tion (0), reduces to unity in the Mestel disk (/3 = 0). 



2.2 Perturbations 

We consider time-dependent first-order Eulerian perturba- 
tions of the equilibria described above and adopt cylindrical 
polar coordinates (R,9,z). Thus, 

E -> E cq (i?) + SZ(R,9,t), 

V> Veq(-R) + dip(R,9,z,t), 

v R -f + 6v R (R,9,t), 

ve — ► v c i TC + 6ve(R,9,t). 

Strictly speaking, ip eq depends upon z as well as R, but we 
have no occasion to evaluate it outside the plane z — 0. 
As the equilibrium is independent of angle 9 and time t, 
the stability analysis is simplified by assuming that the per- 
turbed quantities have the (9, t) dependence characteristic 
of normal modes 



5 oc e 



imQ — iujt 



(13) 



Henceforth, we take this factor as read and often write, for 
example, S~E(R) rather than 5'E(R,9,t). 

The linearized dynamical equations are [cf. Binney & 
Tremaine, 1987, chap. 5] 



-iuSve + ^q$ v R 



mi 

IT 



) 

,5E\ 
E ) ' 



. _<5E r,SvR d . im c 

~ W ~S ~ "r" = -dR 5VR --R 5ve > 
together with Poisson's equation 

V 2 8iP = 4itG5ZS(z). 

We have used the abbreviation 

Co = u> — mQ, 



(14) 



(15) 



(16) 



for the Doppler-shifted frequency in the local rest frame of 
the gas. Taking the curl of the first two of equations (|l4|) and 
using the third to eliminate the derivative of Svr produces 



1 JL 
RdR 



(RSvg) 



-5v[ 



R 



SvrT, 



d 



20.1: 

K 2 



dR \ 2Q.T, 



(17) 



(18) 



The quantity 
- 20E 

is the equilibrium distribution of potential vorticity (or the 
curl of the velocity divided by surface density). Equation 
( p"7j ) states that fluid elements conserve their potential vor- 
ticity in a Lagrangian sense. This is easy to understand phys- 
ically, as the fluid elements must preserve both their mass 
and their circulation (Kelvin's theorem). 

We expect a continuum of modes involving a net change 
in potential vorticity. Such modes will be discontinuous in 
Svg at the corotation radius where Co = 0. This can be any- 
where in the disk. The eigenfrequencies {to} are continuously 
distributed. We are not interested in these modes because 
they can never be unstable (co must be real). If any such 
mode were to grow or decay, the net potential vorticity of the 
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disk would have to change, which is not permitted by our in- 
viscid equations. The issue of interest to us is whether there 
exist growing, non-vortical modes: that is, disturbances with 
the (9, t) dependence ( |l3| ) and vanishing total potential vor- 
ticity. The equations governing the non- vortical modes sim- 
plify in the Mestel disk because £ eq oc R 13 ^ 2 , which is con- 
stant when j3 = so that equation (XL7\) reduces to 



1 8 



(RSve) 



-8vr 



K 02j 
2filT 



0. 



(19) 



RdR K ' R 
Henceforth, we restrict our attention to the Mestel disk. 

3 MODES WITHOUT SELF-GRAVITY 

In this section, on top of the restriction to f3 = 0, we neglect 
the perturbation in the gravitational potential: that is, we 
force Sip — > 0. This is called the "Cowling approximation" 
in the context of stellar oscillations, where it is justified by 
the central concentration of the star's mass or by the short 
wavelength of the modes of interest. Similar justifications 
could be offered for the disk modes in the limit /3 ~ 1, or for 
general /3 in the limit of large m. But we invoke neither limit, 
because we make no pretense of quantitative accuracy until 
Section 4. The true justification for neglecting the gravita- 
tional perturbation is that we thereby obtain equations that 
can be solved exactly. This simplified problem retains the 
elements that make the stability problem ambiguous in the 
self-gravitating case. 

3.1 A Second-Order Differential Equation 

With <W> neglected, the second of equations (|lj) and equa- 
tion (h9() can be solved for Svr and <5£ in terms of the vari- 
able 



y = R5v g 

and its radial derivative: 
dy 



if ioy 

- — mc —-- 

D \ dR 

1 / dy 



Sv r = — 
5E 

E - D { VciTC dR 
where the denominator 

D 



+ mujy 



2 

Weir 



. 2 2 

+ m c 



(20) 



(21) 



(22) 



is constant. Eliminating Svr and 5E from the third of equa- 
tions ([lij) produces a second-order differential equation for 
y, which we write in the equivalent forms 



d^y_ 

dR 2 



+ 



m 
~R 2 



"775" I y = o, 



tfy_ 
dR 2 



(uj - mQ.) 2 - 2fi 2 



rn 



y = 0. 



(23) 



Equation (|23|) is non-singular. In the general (3 ^ case, 
second-order equations can also be obtained for the various 
fluid variables, but at least one coefficient will be singular 
at corotation (il> — > 0) because of the non-zero potential vor- 
ticity gradient. In all other respects, equation (^ ) is typi- 
cal of the general case. The wave is evanescent between the 
Lindblad resonances, Co = ±k or u = mQ ± n. In fact, the 
WKBJ radial wavenumber becomes real somewhat beyond 



the Lindblad resonances (further from corotation) because of 
the lo 2 l<? term. If we had allowed for the self-gravity of the 
mode, propagation would be possible between the Lindblad 
resonances but not in the immediate vicinity of corotation 
unless Q < 1 (see, for example, Binney & Tremaine 1987 or 
Shu 1992) 

Equation ( p3| ) can be reduced to a confluent hypergeo- 
metric function. Let us make the following transformations 
of the independent and dependent variables: 

2iu 



R- 



-R, 

iujR/c 



z/2 



w(z) 



(24) 



where fi is the following dimensionless function of m and the 
Mach number 081): 



i [4(m 2 - 2)M 2 - 4m 2 - l] 



1/2 



(25) 



There should be no confusion between the dimensionless 
variable z and the original cylindrical polar coordinate of 
the same name. Equation (B3h becomes 



Z _ + ( 1 + 2 V . 



> dw . ... 1 
[in + vmM + - 

dz \ 2 



w = 0, (26) 



which is Rummer's equation (Abramowitz & Stegun 1970). 
The solution regular at z — is Rummer's function A/(i + 
i/i + irnM, 1 + 2ifi, z). In fact, if one uses — fi instead of n 
in the transformation (|24|) of the dependent variable, one 
obtains an equation identical to (|3(j) except in the sign of 
fx, whose regular solution is M(| — z/i + imM, 1 — 2i/x, z). 
Hence, the two independent solutions of equation ( pj| ) are 



y±(z) 



a 4±* e -»/a 



± ifj, + imM, 1 ± 2i[i, z 



(27) 



Let us note that Drury (1980) already derived a differential 
equation equivalent to (E3) and used it to find the reflection 
and transmission coefficients for the Mestel disk. He did not, 
however, construct the normal modes, which we now proceed 
to do. 



3.2 Boundary Conditions 

The desired eigenfunction is a linear combination of the in- 
dependent solutions ( p7j ) that satisfies appropriate boundary 
conditions at large and small R. The boundary condition at 
large radius is straightforward. We insist on purely outgo- 
ing disturbances, otherwise it is all too easy to manufacture 
counterfeit instabilities whose exponential "growth" is due 
solely to increasingly powerful transmissions from sources 
at R — oo. One can see directly from the second form of 
equation (^) that the radial wavenumber kR — > ±u/c as 
R — > oo for any value of m provided only that Real(tj) 7^ 0. 
The outgoing wave is proportional to 

+iuiR/c -z/2 

e oc e , 

and the signs of the exponent s are reversed for the ingoin g 
wave. Using standard results ( Abramowitz fc Stegun 1970 ), 
one sees that the asymptotic behavior of y±(z) as \z\ — > 00 
in the cut plane — |7r < arg(z) < |7r is 



y±(*)~r(l±2t/i) 



D «7r( ■jg+imMtihifJ.) ~— im Ai g — z/2 

r(| ±ifi- imM) 
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r(| ± i/x + imM) 



(28) 



We have to cancel the term in e 
combination is 

r(i- «/!)»+(*) r(i + 2i M ) y _( z ) 



so the desired linear 



i^t + imM) r(i + i^i + imM)' 



(29) 



or any fixed multiple thereof. 

The boundary condition at the origin requires more 
thought. First, it is important to realize that an ingoing 
wave reaches the origin - where there is a density cusp - 
in finite time. Propagation to the very centre requires M > 
m/Vm 2 - 2, or equivalently < 9 < (m 2 - 2)/(2m 2 + 2). 
This is never possible for m = 1 disturbances when self- 
gravity is neglected, but can happen in sufficiently cold disks 
for all other azimuthal wavenumbers. The WKBJ dispersion 
relation implied by equation (E3I) is 



(lo — mQ,) 2 



2Q, 2 



in 

c 2 "R 2 ' 

At sufficiently small radius, where Q(R) 2> |w| 
to 



k 



(m 2 - 2)M 2 



2mM , 2, 

-U! + U(U! ). 



(30) 

this reduces 
(31) 



(32) 



R 2 Rc 

The numerator of the first term on the right is almost 4/z 2 [cf. 
eq. (pfj)]. A comparison with the exact solutions ( |2^ ) shows 
indeed that it should be 4/i 2 . Asymptotically, therefore, the 
group velocity is 

duj 2\l 

the upper signs applies if ha > 0, and we have assumed 
m > 2. Thus, the propagation speed is of order the sound 
speed, and an ingoing wave can reach the centre in finite 
time. 

Allowance must therefore be made for a reflected wave. 
By analogy with the situation at large R, we require that 
the origin must neither absorb nor emit wave energy. We 
proceed to implement this principle. 

At sufficiently small R, the Kummer functions M(i ± 
i/i + imM, 1 ± 2i\i, z) — » 1, and we can read off the ingo- 



ing and outgoing parts of the general solution (t 
compatible with the large- R boundary condition: 



that is 



!Jn 



Vout 



r(i + 2ip) 



T(i +ifi + imM) 
T(l - 2i/i) / 2iu 
r(| - ifi + imM) \ c 



R 



(33) 



(In this paper, the symbol "~" means "asymptotically ap- 
proaches" .) Let us define the flux as the total power crossing 
a cylinder centred on the rotation axis. Then, the fluxes car- 
ried by yi n and y out are 



F in 

Cm. 



- Cm 

,(R) 



j(R)\yin\ and F ou 

^iE(i?)Real(oj) 
= M 2 + m 2 ' 



-C m ,Lj(R)\yout 



(34) 



in the limit SI ^> \uj\. The corresponding angular momentum 
fluxes can be obtained by dividing by the pattern speed 
Q p — Real(tj)/m. The coefficient C m ,w{R) can be derived 
by adapting results from Appendix A of Narayan, Goldreich 



-5.5 n 1 — 1 — 1 — 1 — I — 1 — 1 — 1 — 1 — I — 1 — 1 — 1 — 1 — I — 1 — 1 — r 



a 



StfJ 

o 




2 3 4 

Mach number 



Figure 1. Plot of the logarithm of the ratio of growth rate to 
pattern speed s/mf! p as a function of the Mach number M for 
non self-gravitating modes. The diagram shows the dispersion 
relation for modes with azithumal wavenumber m = 2, 3 and 4. 



& Goodman (1987), but for our purposes, all that matters 
is that the ratio of these fluxes is 



fin . 

F ou t |y ou t| 2 ' 

whence we derive the boundary condition 



\yin\- 
\y°ut\ 



1 as R -> 0. 



(35) 



3.3 Dispersion Relation and Growth Rates 



The ratio of the moduli of ingoing and outgoing parts of the 
eigenfunction (133) is 



Vout 



= exp[/i7r — 2(1 arg(tj)] 



Therefore, the inner boundary condition l35h requires 



r (| + imM + in) 



r (| + imM - ifi) 



(36) 



arg(tj) 



tan 

7T 

2 _ 



-r / 



V mQ p 



iln 
4 M 



cosh n(mM + //) 



cosh Tv(mM — /Lt) 



(37) 



Here, the complex eigenfrequency u> has been written as 
mflp + is, where s is the growth rate and O p is the 
patt ern speed of the mode. We have used some identi- 
ties ( Abramowitz & Stegun 1970 ) to reduce moduli of com- 



plex gamma functions to elementary functions in the disper- 
sion relation (^). Fig. 1 shows the variation of the ratio of 
growth rate to pattern speed with the hotness of the disks 
for various azimuthal wavenumbers. For any temperature or 
Mach number, only this ratio is fixed by the dispersion rela- 
tion. The actual values of s or fi p are indeterminate. This is 
because there is insufficient physical information to deduce 
the phase shift of waves reflected from the origin. 
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It may help to give an example of the sort of inner 
boundary condition that would fix the required phase and 
thereby yield a definite discrete spectrum of modes. Suppose 
the origin is enclosed by a rigid wall of small radius Ro <C 
jfcirc/M- Then since no gas penetrates the wall from outside, 
Svr = at R — Ro- Using eq. (£11) to write 5vr in terms of 
y = yin + 2/out, and using eq. (p3) to estimate dy/dR, one 
finds the boundary condition 



2/out (Ro) _ § +mM 



ifi 



yin(Ro) 



+ mM + ifx 



J<t>0 



(38) 



to leading order in Ro <C v c i IC /cu. Notice that the intermedi- 
ate expression above indeed has unit modulus, so that (f>o is 
real: this is a consequence of the fact that a rigid wall absorbs 
no mechanical energy. On the other hand, the asymptotic 
expressions (B3J) yield 



at- ( ^± ] = y.-(,„..VI) H-2/Hu 



2uR 



(39) 



where F is independent of R and u>. Comparing these last 
two equations, we find that 



In \co\ = 



+ 2nn - F(m,M) 
2/i 



ln.Ro, 



(40) 



where n can be any integer. Thus the eigenfrequencies form 
a discrete geometric progression: uJn+i ~ e"'^. The actual 
frequencies depend upon the choice of boundary radius Ro, 
although the ratio of their real to imaginary parts does not. 
For a different kind of physical boundary condition (perhaps 
a free inner edge, for example), one would have a different 
(j>o and hence a different correspondence between Ro and 
the eigenfrequencies. However, arg(u;„) is always given by 
eq. (M), provided uj n < v circ /Ro. 



4 MODES WITH SELF-GRAVITY 

Let us now return to the fully self-gravitating case. In Sec- 
tion 4.1, we derive a recurrence relation between the Mellin 
transforms of the perturbed fluid quantities. This relates 
Mellin transforms with arguments with the same real part 
but differing by integers in their imaginary part. The bound- 
ary conditions for this recurrence relation are discussed with 
some considerable care in Section 4.2 prior to the numerical 
construction in Section 4.3. 



4.1 The Recursion Relation 

As the unperturbed physical quantities vary like powers of 
the cylindrical radius, the Mellin transform (e.g., Carrier, 
Krook & Pearson 1966) is a natural tool for simplifying 
the linear stability equations (jl4|). First, we collect the real- 
space variables into a dimensionless four-component vector 



U(R), (41) 



\ L c c z c 

with the corresponding Mellin transforms u(a) = 
(ui,U2,U3,U4) T defined by 



u(a) 



dR R- 1/2 ~' a U(R). 



(42) 



Both sides have the implicit (6, t) dependence exp(im8 — iujt) 
of course. This is not quite the usual definition of a Mellin 
transform (Carrier, Krook & Pearson 1966), but the factor 
of i? -1 / 2 in the integrand makes for a real value of «o as 
defined by eq. (^6|) below. 

From the second and third of eqs. (|l4|), eq. (|l5|), and 
from eq. (ft9|), we have the following system: 

^-ui(a + i) = —im[Mui(a)+U2(a)] 
ic 

-(ia- ^)u 4 (a), 

—im [A4 112(a) + mi (a) - 

— Mui(a), 
(a) = (M 2 + l)K{a, m)wi(a) 



— U2(a + 1) 
ic 



u 3 (a)} 



ia + -)u2(a) — Mui(a) 



(43) 



Here, K(a, 



Kalnajs 1971; Binney 

ir[|(|- 



is the famous Kalnajs function (Snow 1952; 
Tremaine 1987), 



K(a, m) 



+ m + ia 



)]r[S 



+ m - 



")] 



(44) 



which relates the Mellin transforms of the density and the 
potential perturbations. Note the distinction between Dq = 
27rG£ i?o and u c 2 irc = vf,M 2 /{M 2 + 1). Only the first two of 
eqs. ([I3]) couple different values of the Mellin variable a. So, 
after elimination of U3 and U4 using the last two equations, 
one has a system of the form 



—uia + 1) 

ic 



M(a 

where u = (ui,W2, 

M A 
m 2 



u(a), (45) 
T and Aff (a) is the 2x2 matrix given by: 



M 



a 2 + m 2 + 1/4 



+ M 2 -m 2 (M 2 + 1)K 



MA 
m ^2 



It is helpful to view eq. (|4J) as a recursion relation relating 
Mellin transforms with arguments with the same real part, 
but differing by imaginary increments. We shall refer to such 
an arrangement as a "ladder". 

For neutral modes, u) — » 0, there is no coupling between 
different values of a in eqs. (|43|). Equivalently, there is only 
one rung of the ladder. As realised by a number of inves- 
tigators (e.g., Zang 1976; Lemos & Lynden-Bell 1993; Syer 
& Tremaine 1996; Evans & Read 1998b), exact equiangular 
neutral modes are possible. In fact, a solution exists in the 
form of an exact logarithmic spiral provided det M = 0, or 

[1 - (M 2 + l)K(a, m)][a 2 + m 2 + ~] - M 2 (m 2 - 2) = 0.(46) 

The roots of this equation are real. Let us call them ±«o- If 
we artificially set K — > thereby turning off the self-gravity, 
Qo reduces to the quantity fi defined in eq. (|25[) . 



4.2 Boundary Conditions 

Let us now consider two complementary, semi-infinite lad- 
ders differing in the sign of the real part of the Mellin argu- 
ment: a = +ao + ni and a — — ao + ni, where n is nonneg- 
ative integer that can become arbitrarily large. What are 
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the boundary conditions for the recurrence relation? This 
is a delicate matter and we will consider the bottom and 
top rungs of the ladders in turn. In physical terms, these 
correspond to small- R and large- J? respectively. 

4.2.1 The Top Rung 

At large radii, the WKBJ dispersion relation 



r— J 



2vi 



- V Ak + c 2 k\ 
H 



R ) P 2 #•-.-•-> ( 47 ) 

indicates that self-gravity becomes negligible and the radial 
wavenumber becomes constant, viz., 

(M 2 + l)/2 -mM 



k + 

c 



R 



+ 0(R~ 2 ) if P» ■=2" (48) 



The sign of k has been chosen appropriately for the outgoing 
wave, and the subdominant 0(R~ 1 ) terms have been explic- 
itly retained as they contribute logarithmically to the phase 
of the eigenfunction. It follows from equations ( po| ) - ( |23| ) 
that 

U(R) ~ exp {^>j R H(M 2 +i)/2 ~™M] Uoo (49) 

as R — > 00, where Uoo is another constant vector. This 
can be compared with the eigenfunction in the non self- 
gravitating case [defined by eqs (b3)-(E3)], where the expo- 
nent is only —imM. The term (Mr + l)/2 in eq. (0) stems 
ultimately from the self-gravity in the WKBJ dispersion re- 
lation. At the nth rung of the ladders, the Mellin variable 
has had its imaginary part incremented n times and has 
become 



Q n = ±ao + in 



(50) 



As n — > 00, the Mellin transform (^2|) is dominated by the 
large- R behavior above, so that 



u(a n ) 



w oc r ( s fi 



s n = n + — — i ( ±a + mM 



M' 



(51) 



The direction of the constant vector u-o is determined by the 
relations (po|)-(|2^). The second component of Moo is 0(n~ 1 ) 
smaller than the first component because Svg — y/R ~ 
R~ (<5£/E) asymptotically. As the overall normalization of 
the eigenfunction is arbitrary, the simplest approximation is 



(52) 



In fact, though the WKBJ analysis gives good physical 
insight into what is happening at the top rung of the ladders, 
the boundary conditions ( |5l| ) - jH^ ) are insufficiently accu- 
rate for numerical construction of the modes. This compels 
us to develop an asymptotic analysis to find the higher order 
corrections. The following paragraphs accomplish this task. 
They are primarily matters of detail rather than principle. 
Readers interested mainly in the latter can skip to the next 
sub-section without losing the thread of our argument. 

Rather than deal with two coupled first-order differ- 
ence equations, it proves convenient to recast (EE|) as a 
single second-order difference equation for which standard 
asymptotic techniques are readily available in the literature 
(e.g., Bender & Orszag 1978, chap. 5). Let us write the two- 
component vector u(a n ) as 



u(a n ) 



(Zn,2/n) 



(53) 



Substitution into eq. (^5J) yields a coupled pair of first- 
order difference relations relating (x„+i,y n +i) to (x„,y n ). 
By eliminating the ys, one has a second-order difference 
equation: 



with 

B n 

D n 



A„ + 



B n D n 



B„ 



x„ -\ x n _i = 0. 

B n -i 



(54) 



— (7: + a n - «m ) ~ 0{n) 
rn Z 

a 2 n + m 2 + 1/4 



177! 

m 2 + M 2 - m 2 (M 2 + l)K{a n ,m) 



0(n 2 ) 



M,i 

— (- — a n — «m 
m 2 



tm 

2\ 



0(1) 



0(r 



(55) 



This notation is used because the matrix M n — M(a n ) 
which takes us from the nth to (n + l)th rung of the ladder 
is just 



M n = 



It has determinant 



(56) 



(57) 



For future use, let us also note that the second component 
y n +i can be derived via 

y n+ i = (D n x n+ i — A n x n )/B n (58) 

The coefficients of x n +i,x n and x n -i in (^ij) are O(l), O(n) 
and 0(n 2 ) respectively. Therefore, this suggests the ansatz 
(c.f. Bender & Orszag 1978) 

= Pi a n + Po + P-ia- 1 + P-2a~ 2 + ... (59) 

Xn 

where the {Pk} are independent of n. By substituting this 
ansatz into the difference equation (^) and matching the 
orders of the asymptotic expansion, it is possible to solve 
recursively for the {Pk}- A second-order linear recursion re- 
lation should have two independent solutions, hence there 
are two roots for the leading coefficient: Pi = =pi. Only the 
upper sign is consistent with the required behaviour (^l[), 
and the subsequent terms are then determined. From the 
first two terms, we have 



x„+i 1 
= n H 

Xn 2 



±ao + mM 



M 2 + l 



+ 



+ 



(60) 



Un+i M+im 

Xn+l n 

where the second line has been derived with the aid of (E3) . 
This suggests that as n — > 00, the solution on the nth rung 
approaches 



U(ctn) ~ T(Sn) 



M + imy 
n ) 



(61) 



This is recognised as a more accurate version of eqs ( |5l| ) - 
©• 

From a computational point of view, it is best to carry 
out the calculation (eased by computer algebra) to still 
higher order. We find that: 



© 0000 RAS, MNRAS 000, 000-000 



The Paradox of the Scale-Free Disks 7 



OL n 



(62) 



with 



h = 5 



i mM 



M 2 + l 



(63) 



(2 - n 2 )M 2 - ^ + m 2 + 2i(mX - X 2 - 1) 



+2z/(// + imM) 



(64) 



The second component y n can be found, for example, by 
use of eq. (^j|). These are the boundary conditions that are 
employed in our numerical implementation in Section 4.3 
below. 



4.2.2 The Bottom Rung 

Now, let us see how the bottom rungs of these ladders 
(a = ±ao) correspond to the ingoing and outgoing waves 
far inside the corotation radius. The Mellin transforms have 
poles at a = ±a?o because there is an infinite range of In R 
over which the eigenfunction is asymptotically equal to a 
sum of two power laws. In fact, this representation is exact 
for a neutral mode (uj = 0), where the Mellin transforms 
vanish on the higher rungs of the ladder. In the general case 
uj 7^ 0, one expects that the growing modes are superposi- 
tions of the neutral modes at least for R <C v c i IC /\uj\. In 
other words, 



U(R) ~ U+R 



+ U-R 



-1/2— iao 



(65) 



where the coefficients U± are constant vectors. Hence, if 
the range of integration in (^) is divided into two intervals 
< R < Ro and Ro < R < 00, with Ro <§C Wch-c/M, then 
the contribution from the first interval is approximately 



U+R$ 



-\-ictQ —ia 



U-R~ 



ZOLC\ — %OL 



i(ao — a) i(—ao — a) ' 

provided Imag(a) > 0. Guided by this, we expect that u(a) 
is composed of both a regular part and a singular part with 
simple poles: 



1(a) 



s (±ao) + 



a =F a o 



as a — > ±ao + iO 



(66) 



Here, u± are two-component constant vectors linearly re- 
lated to U±, while u les (a) is a vectorial function that is 
regular at the poles. Furthermore, it must be true that 



M(a )u+ = = M(-a )u~ 
so that u(±Qo + i) is finite: 



(67) 



- i— it(±ao + i) 

c 



lim M (±ao + ie) ■ u(±ao + it) 

dM 

da 



■ u± + M(±a ) ■ w re g(±ao) 
(68) 

Owing to the conditions (|67|), the directions of the vectors 
u± are fixed. Their complex normalization factors are re- 
lated by the requirement that the ingoing and outgoing parts 
of the asymptotic form (pa) have equal moduli. Therefore, 



«+ 



(69) 



In other words, the residues of the poles at a — ±«o are 
related by the boundary conditions at small R. Since the 
ingoing and outgoing wave must carry equal energies, the 
moduli of these two residues must be equal. The relative 
phase of the two residues is arbitrary in the absence of a 
condition on the phases of the ingoing and outgoing waves. 

Thus the large- R boundary condition gives us two com- 
plex constraints, or four real ones, at large n (one complex 
constraint for each of the two ladders). The small- R bound- 
ary condition gives a single real constraint. Our degrees of 
freedom are the complex residues of the poles at the bot- 
tom rung of each ladder (two complex or four real param- 
eters), plus the complex parameter u. Hence we have five 
real constraints and six real free parameters. As in the non 
self-gravitating case, the problem is indeterminate unless we 
impose an additional constraint. This could be the relative 
phase of the residues at iao, or else the real part of ui. 

At the risk of repetition, we now show that just as in the 
non self-gravitating case of Section 3, so here also the ratio 
of moduli of the ingoing and outgoing waves is determined 
by axg(u), while the relative phase is determined by On 
the one hand, we have seen that the complex amplitudes of 
the outgoing and ingoing waves at small 7? are proportional 
to u+ and w_. On the other hand, the recursion relations 
(p4])-(|5^) for the quantities {x n ,y n } defined by eq. ( |H^ ) are 
independent of u>. Therefore the ratio of ingoing to outgoing 
complex amplitudes is 



«4 
U- 



■ n (ao + in) 
t (— Qo + in) 



2ia 



= exp [— 2qq arg(w)] ■ exp [2iao In | 



(70) 



where the unwritten constants of proportionality are inde- 
pendent of lo and n as n — > 00 [cf. eq. (^lj)]. Therefore all 
boundary conditions that neither absorb nor emit energy 
produce the same ratio of real to imaginary parts for all 
eigenfrequencies, but inasmuch as they fix different phase 
shifts at the inner boundary, they produce different discrete 
values for the moduli of the eigenfrequencies. Also, succes- 
sive eigenfrequencies are always separated by n/ao in the 
natural logarithm. 

4.3 Numerical Strategy and Results 

To construct the normal modes, we must adjust the com- 
plex parameter uu so that all the boundary conditions are 
satisfied. For a large value of n, the starting value of u(a n ) 
is fixed on the two ladders using eq. (^i|). The recurrence 
relation ( ji^ ) is successively applied to reduce the imaginary 
part of the Mellin transform variable to unity. This calls 
for the evaluation of the Kalnajs function ( [44] ) for complex 
values of a n with large imaginary parts. The best way to 
proceed is to use Stirling's formula to derive the asymptotic 
result 



K(a n ,m) ~ ± — 
a n 



4m 2 



8ag 



and then use the recurrence relation 



K(a n - 2 ,m) = 



m 2 + (a n - 2 + |i) 2 



K(a n ,m) 



(71) 



(72) 



rn 2 + (a n _2 + ij) 2 

to work down the ladders. The final step to the bottom 
rung warrants further discussion. From eq. (pa), the quan- 
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Figure 2. Plot of the argument of the eigenfrequency (or arctan- 
gent of the ratio of growth rate to pattern speed atan[s/nrf2 p ]) as 
a function of the Mach number JV[ for self-gravitating modes. The 
diagram shows the dispersion relation for modes with azimuthal 
wavenumber m = 1, 2 and 3. 



tity u(±Qo+i) contains contributions from both the singular 
terms (which lie in the null eigenspace) and the non-singular 
terms (which lie in its complement). Let n be the null eigen- 
vector so that M • n = and let p be a projection operator 
that annihilates the non-null eigenvector. Then the residues 
we seek on each ladder are 



u± 



u(±ao 



p T ■ M'(±a ) • n 



(73) 



To satisfy the inner boundary condition, these residues must 
be equal in modulus. The entire numerical algorithm can be 
tested in the non self-gravitating instance by the obvious 
stratagem of setting the Kalnajs gravity factor K(a, m) to 
vanish. Here, the exact answer is already known by virtue 
of the work in Section 3. 

Fig. 2 shows a graph of the ratio of growth rate to 
pattern speed for self-gravitating modes in the Mestel disk. 
Now, growing modes are possible provided the disk is cold 
enough. The condition for instability can be worked out ex- 
actly as 

i 2 + K{0, m)] 



M z > 



-2 + (m 2 + i)A-(0,m) 
where K(0,m) is readily deduced from eq. (^) as 
ir 2 [i(i+m)l 



2 r 2 [i 



+ m 



)] 



(74) 



(75) 



So, for example, m = 1 growing modes are possible only if 
M > 0.869, and m = 2 modes only if M > 0.733. At a fixed 
temperature or Mach number, the ratio of growth rate to 
pattern speed is fixed. But, just as in the non self-gravitating 
case of Section 3.3, the magnitudes of these quantities are 
not fixed as a consequence of the indeterminate nature of 
the phase shift at the origin. 



This paper has resolved the paradox of the scale-free disks. 
By itself, the linear stability problem is ill-posed. Normal 
modes do exist, but the moduli of the eigenfrequencies de- 
pend on the undetermined phase shift with which waves 
are reflected from the centre. For any definite choice of 
this phase, there can be a discrete set of growing modes 
for some disks. On the other hand, the argument of the 
eigenfrequencies — or equivalently the ratio of growth rate 
to pattern speed — is independent of the phase shift. This 
paper has explicitly calculated such ratios with and without 
self-gravity. 

The paradox of the scale-free disks is of interest for 
two reasons. First, as Birkhoff (1960) has indicated, para- 
doxes are of inherent pedagogical value. It sharpens the in- 
sight to track down the flaw in a plausible physical argu- 
ment that nonethless leads to an apparent inconsistency. 
The paradox of the scale-free disks is a paradox of over- 
simplification. In theoretical work, it is customary to develop 
as simple a model as possible to describe objects or phenom- 
ena. The power-law disks, in which all physical quantities 
scale like powers, are attractive candidates for representing 
both galactic and accretion disks. But, they are not suffi- 
ciently rich to act as reasonable models for stability prob- 
lems - unless extra ingredients such as reflecting boundaries 
or cut-outs are added. This is the procedure that was fol- 
lowed in Zang (1976) and Evans & Read (1998a, b) to yield 
normal modes. Second, the paradox draws attention to the 
importance of inner boundary conditions in all models with 
central cusps. The stability properties are crucially affected 
by the manner in which impinging waves are reflected off 
the cusp. Therefore, the correct boundary condition at the 
centre - whether applied directly in a linear stability analy- 
sis or indirectly in a computer simulation - needs careful 
thought. Thinking physically, central cusps are usually a 
consequence of black holes. Probably the most astronomi- 
cally relevant boundary condition is to allow the wave to 
reflect off a central black hole. 
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